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INTRODUCTION 


Transient  solutions  of  the  hyperbolic  type  partial  differential  equation, 
for  example  the  wave  equation  or  the  beam  equation,  are  Important  for  solving 
engineering  problems  such  3S  stress  wave  for  gun  dynamics  or  shock  behavior  of 
penetration  mechanics.  A_t  present  these  equations  are  usually  solved 
numerically  by  the  finite  difference  method  or  by  the  Galerkin  method. 
Considerable  advantage  may  be  obtained  if  the  finite  element  method  can  be 
directly  employed  instead.  Variational  procedures  using  bilinear  formulation 
with  adjoint  variables  can  serve  as  the  theoretical  basis  for  the  derivation 
of  algorithms  using  the  finite  element  method  for  the  hyperbolic  type  partial 
differential  equations  (PDE). 

THE  VARIATIONAL  PRINCIPLE 

A  dynamical  system  can  be  modeled  by  the  following  partial  differential 
equation. 

L(C)  y(C)  =  -Q(C)  (I) 

with  appropriate  boundary  and  initial  conditions.  In  the  above  equation  L  is 
a  linear  operator  in  both  spatial  and  temporal  domain,  y  is  the  dependent 
variable,  Q  is  a  forcing  function,  and  5  represents  all  Independent  variables, 
both  spatial  and  temporal. 

The  inner  product  <  >  of  an  adjoint  forcing  function  Q  and  the  solution 
(y(O)  of  Eq.  (1)  can  be  used  for  the  purpose  of  estimation.  This  inner 


product  is  <Q,y> 


An  accurate  estimation  can  be  made  by  constructing  a  variational 

principle  (ref  1).  By  using  the  adjoint  variable  y  as  a  Lagrange  multiplier 

for  Eq.  (1)  adding  to  <Q,y>,  we  have 
-  A  - 

Jity.y]  -  <Q,y>  +  <y,(Q«-y)>  -  <Q,y>  +  <y,Q>  +  <y,Ly>  (2) 

To  keep  the  system  symmetrical,  let  us  define  the  adjoint  system  as 

L(e)yU)  =  -Q(0  (3) 

By  using  the  original  variable  y  as  a  Lagrange  multiplier  for  Eq .  (3)  adding 
to  <Q,y>,  we  have 

J2ly*y]  *  <Q»y>  +  <y,(^Ly)>  =*  <Q,y>  +  <y,Q>  +  <y,Ly>  (4) 

By  definition,  the  relationship  of  the  adjoint  system  to  the  original  system 
is 

A  - 

D  =■  <y,Ly>  -  <y,Ly>  =*  0  (5) 

where  D  Is  the  bilinear  concomitant  (ref  !)•  Combining  Eqs.  (2),  (4),  and  (5) 


one  obtains 


*  J2 


In  order  to  keep  the  functional  symmetrical,  we  have 

A  1 

J  "  "  lJl  +  J2l  < 

which  is  of  the  form 

1  -  1  — 

J  -  <Q,y>  +  <y,Q>  +  -  <y,Ly>  +  -  <y,Ly>  ( 

2  2 

To  show  that  the  above  functional  satisfies  both  the  original  and  the 


adjoint  systems,  let  us  take  the  first  variations  of  Eqs.  (5)  and  (6)  which 
gives 


6J  »  6j(<Sy)  +  <5j(  6y) 


(7a) 


^Stacey,  Weston,  M.  Jr.,  Variational  Methods  in  Nuclear  Reactor  Physics, 
Academic  Press,  1974.  *  . " 
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where 


6J(6y) 

1  -  1  -  - 
■  <6y,Q>  +  -  <6y,Ly>  +  -  <y,L5y>  •  0 

(7b) 

and 

5j(3y) 

»  <<5y,Q>  +  -  <5y,Ly>  +  -  <y,L6y>  =  0 

2  2 

(7c) 

Also 

where 

6D  »  6D(  6y)  +  6D(  6y) 

(3a) 

and 

6D(6y)  -  <5y,Ly>  -  <y,L<5y>  *■  0 

(8b) 

5D(6y)  “  -  <6y,Ly>  +  <y,L<5y>  ■  0 

(8c) 

From  Eqs. 

(7b)  and  (8b) 

we  obtained 

1  -  1 

6j(6y)  -  <6y,Q)  +  -  <<$y,Ly>  +  -  <6y,Ly>  -  <6y,(Q+Ly)>  -  0 

(9) 

For  arbitrary  <5y  satisfying  certain  general  limitations  on  the  boundaries  It 
can  be  shown  that  the  Euler-Lagrange  Equation  for  the  original  system  In  Eq. 
(1)  is  satisfied.  From  Eqs.  (7c)  and  (8c)  we  get 


6j(6y)  -  <6y,Q>  -  -  <5y,Ly>  +  -  <«y,Ly>  -  (6y,(Q+Ly)>  -  0  (10) 

2  2 

For  arbitrary  variation  6y,  the  Euler-Lagrange  Equation  for  the  adjoint  system 
In  Eq.  (3)  is  also  satisfied. 


INTEGRAL  OF  BILINEAR  EXPRESSION 

The  Integral  of  a  bilinear  expression  for  a  two-dimensional  problem 
having  second  order  partial  derivatives  In  time  and  fourth  order  partial 
derivatives  in  space  can  be  written  as 

I  *  /  /  £J[y(x,t)  ,y(x,t)]dtdx  (11) 

*o  to 


3 


where  ft[y,y]  is  a  given  bilinear  expression  in  the  form 

fl(y.y)  =■  ytyt  _  w2yy  -  a2yxyx  ~  b2yxxyxx  (12 

The  subscripts  t  and  x  Indicate  the  partial  derivatives  for  the  functions  y 
and  y. 

Equation  (12)  can  be  Integrated  by  parts.  Two  different  forms  of 
Integration  and  end  conditions  are  given.  The  first  form  of  the  Integral  Is 
obtained  by  integrating  by  parts  on  the  adjoint  variable. 


tb  xb  _  xb  _  tb 

X1  m  -J  I  yLydtdx  +  J  yty|  dx 


to  *o 


*o  uo 


tb  -  Xb  _  Xb  -  Xb 

+  J  {-b2yXxyxl  +  (b2yxx)xy|  -  a2yxy|  >dt 


(13a) 


On  the  other  hand,  we  can  perform  integration  on  the  original  variable  to  give 


tb  xb  __  xb  _  tb 

12  -  “J  J  yLydtdx  +  J  yty|  dx 
tQ  xQ  x0  tQ 


tb  -  xb  -  xb  xb 

+  J  t-b2yxxyxl  +  b2(yxx)xy|  -  a2yxy|  >dt 


(13b) 


To  keep  the  form  symmetrical,  we  take  the  average  of  the  above  two  expressions 


1  1  ,xb  bb  l  -  —  1  xb  tb 

1  “  T  X1  +  T  *  ~J  J  “(yLy+yLy)dtdx  +  -  J  (yty+yty)l  dx 


xo  tQ  2 


1  tb  -  -  xb 

+  -  /  (-a2)(yxy+yxy) |  dt 


1  tb  -  -  xb  i  tb  -  xb 

+  -  /  (-b2)(yxxyx+yxxyx) I  dt  -  -  /  [(-b2yxx)xy  +  (-b2y*x)xy]  dt  (i4) 


Ly  -  ytt  +  w  y  -  a*yxx  +  b*yxxxx 


(15a 


and  _  .  _ 

Ly  -  ytt  +  uty  -  a2yxx  +  b2yxxxx  (15b! 

For  a  fourth  order  spatial  partial  and  a  second  order  temporal  partial  system 

Eq.  (5)  becomes 


*b  «b  -  *b  ,cb  — 

D  ■  J  J  yLydtdx  -  J  J  yLydtdx 


(16a; 


xo  co  *o  co 

By  equating  Eqs.  (13a)  and  (13b)  and  solving  for  D  in  Eq.  (16a)  we  are 
converting  the  double  integral  into  single  Integrals  in  terms  of  the  boundary 
conditions* 

We  can  express  the  quantity  D  as  the  sum  of  three  parts  for  end 
conditions  Dj_t  I>2,  and  D3  as 


D  -  Dx  +  D2  +  D3  (16b) 

The  terms  in  Involve  the  Initial  conditions  of  y  and  y  as 

*b  -  tb  -  tb 

Di  *  /  <yty|  -  7t/l  >dx 

xo  to  co 


Di  -  /  dx{[yt(x,tb)y(x,tb)  -  yt(*»tb)y(x,tb)] 
xo 

-  [ytU.toJyC^.to)  “  7t(x»to)y(x»to)l}  (17a) 

The  terms  in  D2  involve  the  boundary  conditions  from  the  second  partlals  of  y 
and  y  as 

bb  —  xb  —  xb 

D2  -  J  (-a2)(yxy|  -  y^l  >dt 

t0  Xo  Xq 

,fcb  - 

D2  -  J  dt {-a^[yx(xb, t)y(xbl t)  -  yx(xb,t)y(xb,t) 
co 

+  a2tyx(*o»t)y<xo»t)  -  7x(xo»t)y(x0,t)] } 


(17b; 


The  terms  in  D3  involve  the  boundary  conditions  from  the  fourth  partials  of  y 

and  y  as 

D3  -  /  b  {-b2yxxyx|  b  +  b2yxxyxJ  b  +  (b^xxJx/l  b  +  <-b2yxx)xy|  b}dt 


D3  -  /  dt{-b2[yxx(xb,t)yx(xb,t)  -  yjo^xh.Oy^Xh.t)  ] 
co 

"*■  b2[yxx(x0,t)yx(x0,t)  -  yxx(x0»t)yx(x0,t) ] } 


+  /  dt(-b2[-yxxx(xb,t)y(xb,t)  +  yXxx(xb»t)y(xb,t)] 
fco 

+  b2[-yxxx(x0,t)y(x0,t)  +  yx(x0,t)y(x0,t)] } 

In  order  that  D  =  0  in  Eq.  (16b)  it  is  sufficient  that 


(17c) 


=  0 
D2  =  0 
03  =  0 


(18a) 

(18b) 

(18c) 


THE  SYMMETRICAL  ADJOINT  SYSTEM 

The  adjoint  independent  variable  T  in  Figure  1  can  be  expressed  as 

Tb  -  T  t  -  t0 


which  gives 


Tb  -  To  tb  -  t£ 


x  -  xb  for  t  ■  tc 


-  t0  for  t  -  tb 


(20a) 

(20b) 


It  is  noted  from  Eq.  (19)  that 


Tb  “  To  -  tb  -  tc 
X  -  Tb  +  t0  -  t 


(21a) 

(21b) 


dx  ■  -dt 


(21c) 
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(  2  Id ) 


d  d 

dt  dt 


and 


y(x,t)  =*  y(x,T  =  xb+t0-t) 

(21e) 

Let  us  assume  that 

the  adjoint  system  shown  in  Figure  l  gives 

y(x,t=t)  *  y(x,t=tb+t0-t) 

(22a) 

yt(x,t«t)  =  -yt(x,t=tb+t0-t) 

(22b) 

yx(x,t=t)  *  yx(x,t=tb+t0-t) 

(22c) 

where  t  is  a  dummy 

variable  for  t. 

We  may  define  the  adjoint  system  as  the  Image  reflection  In  the  time 
domain  of  the  original  system.  Equation  (22)  yields  the  following  known 
initial  conditions 

y(x,t=tb)  =  y(x,t=t0)  (known)  (23a) 

A  A 

yt(x,t«tb)  =  ~yt(*>t=to)  (known)  (23b) 

The  interpretation  of  the  above  equations  gives  the  Initial  conditions  of  the 
original  system  as  the  far  end  conditions  for  the  adjoint  system,  since  the 
adjoint  system  is  a  reflected  mirror  of  the  original  system  in  time. 

INITIAL  CONDITIONS  FOR  THE  ADJOINT  SYSTEM 

We  take  a  symmetry  approach  for  the  initial  conditions  of  the  adjoint 
system  as 


y(x,t=*tb) 

“  y(x, t=t0) 

,  yt(x,t»tb) 

=  -yt(x,t*t0) 

(24) 

y(x,t-t0) 

*  y(x,t=tb) 

,  yt(x,t-t0) 

a  -yt(x,t=tb) 

(25) 

Thus  Eq.  (17a)  becomes 
*b 

°1  m  I  dx{[yt(x,t»tb)y(x,t=*t0)  +  yt(x,  t=tQ)y(x,  t-tb) 

*o 

-  fyt(x,t«t0)y(x,t»tb)  +  ytCx,t=*tb)y(x,t=t0)]  }  =  0  (26) 
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Since  the  Integrand  of  Eq.  (26)  is  zero,  the  above  satisfies  Eq.  (18a).  The 
initial  conditions  in  Eq.  (25)  are  given.  Therefore 

6y(x,t-tij)  =*  6y(x,t=t0)  -  0  (27a) 

<Syt(x,t-tb)  "  -6yt(x,t-to)  “  0  (27b) 


THE  GENERALIZED  BOUNDARY  CONDITIONS 

Let  us  consider  the  operator  L  in  Eq.  (15a)  for  two  different  cases  as 
follows. 

A.  For  the  wave  equation  we  have 

Ly  *  ytt  -  a2yxx  (28) 

Let  us  assume  that  elastic  springs  are  installed  at  the  ends  such  that 

YxUb.t)  “  kbY(xb,t)  *  yx(xb»c>  *  XbYCxb.t)  (29a) 

YxUo.t)  =  -k0y(x0,t)  ,  yx(x0,t)  “  ~lt0y(x0,t)  (29b) 

This  is  equivalent  to  state  that  the  fixed  end  condition  for  a  prismatic  bar 
is  kb  -  k0  00  and  the  free  end  condition  is  *  k0  +  0.  If  Eq.  (29)  is 
substituted  into  Eq.  (17b)  we  have 

D2  ■  0 

B.  For  the  beam  equation  we  have 

Ly  =  ytt  +  b2yxxxx 
Two  sets  of  springs  are  Incorporated  at  the  ends.  They  are: 

(1)  Torsional  springs  relate  the  moments  (the  second  partials)  with  the 
slopes  (the  first  partials) 

Yxx^b^)  ■  %yx(xb*t)  yXx(xb>t)  ■  ^x^h^)  (32a) 

yXx(x0,t)  -  -n0yx(x0,t)  yXxCx0,t)  -  -n0yx(x0,t)  (32b) 


(30) 


(31) 
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(2)  Linear  springs  relate  the  shears  (the  third  partLals)  with  the 


deflection  (no  partials) 

yxxx^b.t)  "  Cby(xb,t;  yXXx(*b>c)  “  cb/<xb.O  (33a) 

yxxx(xo.O  -  -C0y(x0.c)  yxxx(xo.t)  -  -C0y(x0,t)  (33b) 

By  substituting  Eqs.  (32)  and  (33)  into  Eq.  (17c)  we  have 

D3  =  0  (34) 

Table  I  shows  the  assignment  of  the  spring  constants  for  various  physical  end 
conditions. 


TABLE  I.  GENERALIZED  BOUNDARY  CONDITIONS 


At  Fixed  End 

At  Hinged  End 

At  Guided  End 

At  Free  End 

y**y=0 

y=y=0 

yxayxa° 

yxxayxxa° 

yx=yx=° 

yxx=yxxa° 

yxxx=yxxxa° 

yxxxayxxx=° 

6y=6y=*0 

6y*6y=*0 

6yx=6/xa° 

6yxxa6yxxa° 

Syx^yx3*0 

6yxxa6yxxa0 

6yxxxa<5yxxxa0 

6yxxxa6yxxxa° 

Torsional  Spring 

A 

Yxx^Yx  n  +  ” 

n  ♦  o 

n  00 

n  -*•  0 

Deflection  Spring 

yxxxacy  c  >  °° 

c  -►  to 

c  ♦  0 

c  +  0 

Spring 

yx=*ky  6y»6y=0 

k  00 

k  ♦  0 

k  *  undetermined 

6yx“6yxa° 

THE  FIRST  VARIATION 


where 


The  sum  of  the  two  functionals  Is  obtained  by  adding  Eqs.  (6c)  and  (14) 


*b  tb  -  - 

J  +  I  -  J  J  (Qy+yQ)dxdt  +  T  +  W  +  B 
xo  fco 

1  xb  -  -  tb  1  ,tb  -  -  xb 

t  -  -  J  (yty+yty)l  dx  »  w  ■  -  J  (-a2)(yxy+yxy)l  dt 

Z  Xq  tQ  2  tQ  Xq 


1  ^b  ,  -  -  ,*b  1  tb  ,  -  --  ,*b 

B  -  r  J  (“b  )<yxxyx+yxxyx)l  dt  -  -  /  [(-b2yxx)xy  +  (-b2yxx)xy|  dt  (36) 


By  taking  the  variations  6y  and  Sy  separately,  we  let 

SJ  -  6 J ( 6y)  +  6j(6y) 

Then  one  obtains  from  Eqs.  (35)  and  (36)  that 


$J(6y)  «  -6l(6y)  +  //  Q6y  dxdt  +  6l(5y)  +  6W( 6y)  +  6B( 6y)  *  0 


where 


-  1  Xb  -  ,tb  1  tb  -  ,xb 

«T(6y)  -  -  /  (yt^y+ydyt)!  d*  »  5W(6y)  -  -  /  (-a2)(yx6y+y6yx)|  dt 


2  xo 


6B(6y)  -  \  J  b  (“b2)(yxx6yx+yx6yxx)|  bdt 


-  -  /  b  [(-b2)yxxx6y  +  (-b2)ySyxxx]  b  dt 


where  -6l(5y)  can  be  derived  from  Eqs.  (11)  and  (12)  as 
-  xb  tb 

-6l(6y)  -  -/  /  (yt6yt"“2yiSy”a2yx6yx"b2yxx6yxx)dxdt 

xo  t0 

The  second  term  on  the  right  side  of  Eq.  (37)  is 

5j(5y)  -  -6l(6y)  +  / /  Q6y  dxdt  +  ST(5y)  +  6W(5y)  +  6B(Sy)  -  0 


where 


1  ,xb  - 


cb 


l  ,cb 


,  xb 


6T(6y)  -  -  /  (yt^y+ySyt)  I  dx  ,  5W(<5y)  =  -  /  (-a2)(yx6y4-y6yx)  |  dt 


2  x, 


2  tr 


1  tb 


,  xb 


«B(6y)  -  -  /  (-b2)(yxx6yx+yx6yxx)|  dt 
2  to  xo 


and 


t(-b2yxxx)<sy-b2y<Syxxx] 


(40) 


xb  tb  -  - 

-6i(6y>  -  -  /  /  (yt6yt-a)2y(Sy“a2yxi5yx(“b2)yxx6yxx)dxdc  (41> 

xo  ho 

It  is  noted  that  Eqs.  (38)  and  (40)  are  exactly  the  same  form,  where  Eqs.  (39) 
and  (41)  are  also  similar. 

For  the  beam  equation  it  is  noted  that  the  high  partials  in  Eqs.  (38)  and 
(39)  can  be  replaced  by  Eqs.  (32)  and  (33).  The  variations  of  the  adjoint 
higher  partials  from  these  equations  can  be  written  as 

Syxx^b.t)  “  Myx^bjO  5yxxx(xb>h)  **  cb6y(xb,t)  (42a) 

6yXx(xo.t)  -  -n0dyx(xo.t)  6yxxx(xo.t)  =  -c06y(x0,t)  (42b) 

Equations  (38)  and  (39),  with  the  aid  of  Eqs.  (32),  (33),  and  (42),  are 
the  key  equations  to  be  used  for  the  finite  element  method.  It  is  noted  that 
the  first  variation  6j(6y)  is  the  same  as  the  first  variation  6J( 6y)  by  adding 
or  dropping  the  bar  on  top  of  the  variables  and  their  variations.  We  do  not 
need  to  solve  for  the  adjoint  system  in  Eqs.  (40)  and  (41)  since  they  give 
exactly  the  same  solutions  as  that  of  the  original  system. 
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SECOND  VARIATIONS 


V 


xb 

dt 

(45c) 

xo 

xb 

;)l  dt 

(45d) 

62B  -  -  /  b  (-b2)(6yxx6yx+6yx6yxx)|  dt 


1  tb  -  -  xb 

+  ~  I  b2(<Syxxx6jM-6y6yxxx)  |  dt 

2  tQ  xo 

1  cb  o  -  xb 

62W  »  -  J  (-a2)(6yx6y+6y6yx) |  dt 


The  second  variation  of  I  is  obtained  from  Eqs.  (39)  and  (41)  as 

?  U2  a2i 
S2T  u  (---) 

2  u=o 

=  -  6y[  6l(<5y)  ]  +  -  6y[6l(«y)] 

2  2 

1  ,xb  tb  ,  -  ,  , 

-  -  /  /  (Syt6yt_w  6y6y~a  5yx6yx-b  6yxx6yxx)dxdt 

2  xQ  tQ 

1  xb  tb 

*■  -  f  I  (5yt6yt-u26y6y“a26yx,syx-b26yxx6yxx)dxdt 

2  xQ  t0 

xb  tb 

62i  -  /  /  (6yt«yt"w26y5y-a2(Syx6yx~b26yxx6yxx)dxdt  (45e 

xo  co 

Substituting  Eq.  (27)  into  Eq.  (45b)  we  have 

62T  -  0  (46a 

For  all  the  end  conditions  in  Table  I  either  the  variations  Syxx  and  6yxx  mus 
vanish  or  Syx  and  6yx  must  vanish.  Thus,  the  first  term  on  the  right  side  of 
Eq.  (45c)  is  zero.  Similarly,  for  all  the  end  conditions  in  Table  I  either 
the  variations  6yxx  and  6yxx  must  vanish  or  6y  and  5y  must  vanish.  Thus  the 
second  term  on  the  right  side  of  Eq.  (45c)  is  also  zero.  The  third  term  is 
zero  except  at  the  guided  end.  Thus,  In  general 


^  *  o 


In  Table  I,  except  the  free  end,  either  the  6yx  and  6yx  must  vanish  or  6y 
and  6y  must  vanish.  Thus,  one  obtains 

62W  ■  0  (46c) 

This  reduces  the  second  variations  <$2J  to 

62J  -  -S2J  (47) 

as  given  in  Eq.  (45e). 

Substituting  Eq.  (45e)  Into  Eq.  (47)  gives 
tb  xb 

<52J  »  /  /  [-5yt(x,t)6yc(x,t)  +  oo26y(x, t) 6y(x, t)  + 

to  xo 

+  a26yx(x,t)6yx(x,t)  +  b26yxx(x,t) 6yxx(x,t)]dxdt  (48) 

In  order  that  the  functional  J  is  an  extremum  (refs  3,4),  the  second  variation 
62J  must  be  either  positive  semi-definite  or  negative  semi-definite,  l.e., 

62J  >  0  (or  62J  <  0)  (49) 

The  above  is  a  necessary  condition  for  a  minimum  (or  a  maximum). 

The  adjoint  variations  in  Eq.  (48)  may  be  obtained  by  the  relations  given 


in  Eq.  (22)  as 

a  ^ 

6y(x,t-t)  -  6y(x,t-tb+t0-t)  (50a) 

$yt(x>t”t)  -  -6yt(x,t-tb+t0-t)  (50b) 

SyxC^t-t)  -  6yx(x,t-tb+t0-t)  (50c) 

The  variations  of  adjoint  initial  conditions  can  be  derived  from  Eq.  (23)  as 

A  A 

<5y(x,t"tb)  “  6y(x,t*tb)  ■  0  for  all  x  (5la) 

6yt(x,t-tb)  -  -$yt(x,t-t0)  =*  0  for  all  x  (5lb) 


^Gelfand,  I.  M.  and  Formin,  S.  V.,  Calculus  of  Variations,  Prentice-Hall, 
1963. 

^Sagan,  Hans,  Introduction  to  the  Calculus  of  Variations,  McGraw-Hill,  1969. 
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By  substituting  Eq.  (51)  into  Eq.  (48),  we  have 


where 


62J 


*b 

j  P(x,t)dxdt 


P(x,t)  =  6yt(x,t)6yt(xttb+t0-t)  +  u)26y(x ,  t )  6yx(x ,  tb+tQ- 1 ) 
+  a2<5yx(x,t)6yx(x,tb+t0-t)  +  b26yxx(x,t)  6yxx(x,tb+t0-t) 


(52a) 


(52b) 


SENSITIVITY  RELATIONSHIP 

In  order  to  show  that  the  second  variation  of  the  functional  .1  Is 
positive  seml-definlte,  one  needs  to  obtain  the  variations  of  the  function  and 
Its  partlals  together  with  that  of  the  adjoint  functions  and  its  partlals  as 
indicated  In  Eq.  (48).  We  can  get  these  variations  through  the  study  of  the 
sensitivity  coefficients  (ref  5)  and  Its  relationship  to  the  parameters  given 
in  Eq.  (43).  Let  the  forcing  function  In  Eq .  (1)  be 

Q(x,t)  -  qf(x,t)  (53) 

It  Is  assumed  that  the  forcing  function  parameter  q  is  subject  to  a  small 
constant  perturbation  <5q  as 

q  =*  qQ  +  6q  (54) 

Then  the  variation  of  the  function  y  is 

3y(x, t) 

6y(x,t)  - - 6q  =»  v(x,t)6q  (55a) 

3q 

where 

3y 

V(x,t)  -  “  (55b) 

3q 

The  quantity  v  is  the  sensitivity  coefficient  for  the  variation  Sy(x,t)  due  to 
a  small  constant  perturbation  6q. 

^Toraovlc,  Rajko,  Sensitivity  Analysis  of  Dynamic  Systems,  McGraw-Hill,  1963. 


The  original  PDE  in  Eq.  (15a)  can  be  written  aa 


ns 


6  *  Ly  +  Q 

*  ytt  +  “2y  “  a2yxx  +  b2yxxxx  +  qf(x,t)  =*  0  (56) 

Due  to  the  perturbation  of  q  the  change  of  <|>  obeys  the  following  relationship 

3 6  3ytt  36  9yxx  36  3yxxxx 

- + - + - +•  f(x,t)  -  0  (57) 

3ytt  3<i  3yxx  3<i  3yxxx  3fi 

It  ta  also  noted  from  Eq.  (56)  that 


36  36  9 

- a  1  > - =  U>* 

3ytt  3y 

36  2  36  2 

3yxx  3yxxxx 

Using  the  definition  in  Eq.  (55b)  the  partials  can  be  Interchanged  as 

3ytt  32  3y 

3q  3  3t2  (3q^  =  Vtt 

3yXX  32  3y 

~3q  “  3x2  ( 3q)  *  Vtt 


(58a) 


(58b) 


(59a) 


(59b) 


3yXXXX  3 4  3y 
"Iq"  =  3?  (3q}  “  Vxxxx 


(59c) 


Substituting  Eqs.  (58)  and  (59)  into  Eq.  (57)  we  have 

^tt  ”  ®2^xx  b2vxxxx  ^  ^(x*b)  *  0  (60) 

If  we  compare  the  definitions  of  variation  in  Eq.  (43a)  with  the  definition  of 
sensitivity  relationship  in  Eq.  (55a)  we  have 

<5y(x,t)  -  wn(x,t)  -  (<$q) v(x,t)  (61) 


which  gives 


u(x,t)  -  v(x,t) 


(62a) 


and 


<$q  =  P 


(62b) 


Thus  Eq.  (60)  becomes 

ntt  +  U)2n  -  a2nxx  +  b^rixxxx  +  f(x,t)  =*  0  (63) 

which  gives  the  PDE  of  the  variations  of  the  original  system. 

If  we  compare  Eq.  (63)  with  Eq.  (56)  we  see  that  the  variation  n(x,t)  ** 
p-16y(x,t)  in  Eq.  (63)  takes  the  place  of  the  function  y  in  Eq.  (56)  with  q  = 
l.  Therefore,  the  POE  for  the  variations  is  unchanged  except  by  a  scale 
factor.  Thus  the  solution  of  the  variation  6y(x,t)  has  the  same  form  as  that 
of  the  original  function  y. 

Similarly  for  the  adjoint  system  one  can  obtain 

<5y(x, t)  =»  pn(x,t)  -  (Sq)v(x,t)  (64) 

n(x,t)  =  v(x,t)  (65a) 

6q  =*  p  (65b) 

and 

ntt  +  w2n  -  a2nxx  +  b2nxXXx  +  f(x»t)  “  0  (66) 

which  is  the  PDE  of  the  variations  of  the  adjoint  system. 

EXTREMAL  OF  FUNCTIONAL  FOR  A  SIMPLE  OSCILLATOR 

To  show  that  <$2J  must  be  positive  semi-definite  we  start  with  an  example 
by  a  simple  harmonic  oscillator  with  no  forcing  function.  Thus  from  Eq.  (63) 
we  have  the  ordinary  differential  equation  (ref  6) 

ntt  +  uj2n  =  0  (67) 

^Shen,  C.  N.  and  Wu,  Julian  J.,  "A  New  Variational  Method  for  Initial  Value 
Problems,  Using  Piecewise  Hermite  Polynomial  Spline  Functions,"  ARO  Report 
81-3,  Proceedings  of  the  1981  Array  Numerical  Analysis  and  Computers 
Conference,  1981. 
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The  solution  for  the  above  equation  Is 


Sy  ■  un  3  A  cos({Jt+9) 
6yt  ■  *  -uA  8ln(u)t+6) 


(68a) 


(68b) 


Both  A  and  0  can  be  determined  from  the  following  given  initial  conditions 


6y(t»0)  «  6y(0)  ■  A  cos  9 
6yt(t*0)  ”  6yt(0)  “  -wA  sin  6 


(69a) 

(69b) 


For  computation  by  the  finite  element  method  the  increment  time  Is  taken 


as  T  which  gives 


n  n 

T  -  tb  -  tc  -  (-)(-) 


where  n  *  1,2,3. . . 

The  image  function  becomes 


6y(t*T-t)  **  A  cos[  0+u>(T-t)  ] 

a 

6yt(t*T-t)  *  -wA  sinf  0+tj(T-t)  ] 


(71a) 

(71b) 


For  the  ordinary  differential  equation  we  have  the  second  variation  from  Eq. 
(52)  which  gives 

T 

62J  -  /  [6yt(x,t»t)6yt(x,t*T-t) 

o 

A  A 

+  w2Sy(x,t“t) 6y(x,t*T-t) ]dt  (/ 

Separating  Eq.  (72)  into  two  parts  and  using  Eqs.  (68)  and  (71)  we  have 

62J  -  62J[5yt]  +  <52J[w5y]  (73a) 


where 


UJi 

62J[6yt]  ■  /  u2A2  sin(  9+ut)8in(  0+coT-o)t)d(  <*>t) 


(73b) 


62J[u>6  ]  -  /  w2A2  cos(  9+u>t)cos(  9+ufT-u>t)d(  wt) 


(73c) 
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k  ■  i  iin  i  i  i 


oTT  =*  n(-) 
2 


which  is  a  multiple  of  it/2. 


(73d) 


The  trigonometric  relationship  for  Eq .  (73)  is 

sln(wt+3)  =  sin  wt  cos  6  +  cos  wt  sin  9 
cos(wt+9)  =■  cos  wt  cos  9  -  sin  wt  sin  9 
sin( 9+wT-wt)  *  ~sin[wt  -  ((H-nw/2)] 

»  -sin  wt  cos(  (H-nir/2)  +  cos  cut  sin(  9+tiit/2) 
and  cos( 9+wT-wt)  *  cos[wt  -  (  O-t-mi/2 )  ] 

*  cos  wt  cos(9+nir/2)  +  sin  wt  sln(9+nn/2) 
For  the  case  when  n  is  odd,  we  have 


cos(  9+nti/2)  *  (-1)  2  sin  9 
n-1 

sin(  &fnir/2)  >•  (-1)  2  cos  q 
For  the  case  when  n  is  even,  we  have 

cos(9+nn/2)  *»  (-l)n/2  cos  9 
sin(9+nn/2)  »  (-l)n/2  sln  6 

First,  we  take  the  case  when  n  is  odd.  Substituting  Eqs.  (74)  and  (75) 
Eq.  (73),  one  obtains 

62J[<5yt]  »  w2A2/  {(sin  wt  cos  9  +  cos  wt  sin  9) 
o 

•  [-  sin  u>t(-l)(n‘*‘l)/2  sin  9  +  cos  wt(-l)(n~^)/2  cos  9]}d(wt) 


(74a) 

(74b) 


(74c) 


(74d) 


(75a) 


(75b) 


(76a) 

(76b) 


$ 

s 

<»■ 

I 

$ 

* \ 


(_l)(n~l)/2  u2A2/  [sin  9  cos  9  +  sin  wt  cos  wt]d(wt) 
o 

■  (-l)(n-l)/2  w2A2[-  +  —  sin  9  cos  9] 

2  2 


(77a) 


*  *  *  *  -  *  »  *  «  *  ^  t'k  iN  •  A  % 


6zJ[u)5y]  »  u)zAzJ  {(cos  u)t  cos  0  -  sin  <ot  sin  0) 
o 

[cos  sin  0  +  sin  tot  (-l)(n-l)/2  cos  0]}d(wt) 


-  (-DCn-D/Z  W2A2 {.sin  0  cos  0  +  sin  cot  cos  uot]d(u>t) 

o 

=  (-l)(n-l)/2  o)2A2  [-  -  — )  sin  0  cos  0] 

2  2 

From  Eq.  (73a)  when  n  is  odd  we  have 

52  J  -  (-i)(n-l)/2  u)2A2 {[-  +  --  sin  0  cos  0]  +  [-  -  —  sin  6  cos  0] } 

2  2  2  2 


(77b) 


62J  =  ( —  1 )  ( 1 ) / 2  (i)2A2 


(77c) 


In  particular  for  n  »  1,  one  obtains 


fi2J  =  u>2A2  >  0 


(78a) 


which  gives  a  minimum  for  the  functional  J.  For  n  =*  3 

S2J  -  -w2A2  <  0  (78b) 

which  gives  a  maximum  for  the  functional  J.  It  is  noted  that  S2J  is 

Independent  of  0  which  is  related  to  the  starting  conditions.  It  is  also 

independent  of  the  polarity  of  A  since  it  appears  in  terms  of  A2. 

Now  we  take  the  case  when  n  is  even.  Substituting  Eqs.  (74)  and  (76) 

into  Eq.  (73),  one  obtains 

,  „  nn/2 

S*J[6yt]  ■  u>6A {(sin  o>t  cos  0  +  cos  wt  sin  0) 
o 

[-  sin  ut  (-l)n/2  C03  0  +  cos  ut(-l)n/2  gin  0]}d(u)t) 

■  (-l)n/2  u2A2  /  [-  sin2  u>t  cos2  0  +  cos  2  u>t  sin2  0]d(u>t) 


(_l)ii/2  oj2a2(-  cos2  0  +  sin2  0)mr/4 


(79a) 
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62J[u)5y]  =  uj2A2/  {(cos  u)t  cos  0  -  sin  wt  sin  0) 
o 

•  [cos  u>t(-l)n/^  cos  9  +  sin  ut(-l)n/2  sin  9]}d(wt) 

=  (-l)n/2  oi2A2  /  [cos2  wt  cos2  6  -  sin2  u>t  sin2  6]d(wt) 
o 

=  (-l)n/2  w2A2(cos2  8  -  sin2  0)nTr/4  (79b 

From  Eq.  (73a)  when  n  is  even  we  have 

52J  =  (-l)n/2  id2A2{(-  cos2  8  +  sin2  8)  +  (cos2  6  -  sin2  9)  }nit/4 

S2J  =■  0  for  all  n  =  even  (79c 

We  can  conclude  here  that  the  functional  J  definitely  (ref  6)  has  a 
minimum  if  wT  ®  x/2,  or  T  is  a  quarter  of  the  natural  period  of  the 
oscillation  t  =*  2tt/uj .  Moreover,  from  Eq.  (70)  for  n  =  1 

T  =  cb  “  co  a  ^/C2^)  51  T/4  (30a 

If  n  =  2  and  <52 J  =  0  in  Eq.  (79c),  we  have 

T  =*  tjj  -  t0  <  n/ uj  =»  t/2  (80b 

This  is  the  upper  limit  of  the  increment  we  chose  for  T,  above  which  the 
minimum  of  the  functional  J  is  not  guaranteed. 

EXTREMAL  FOR  A  SIMPLY- SUPPORTED  BEAM  WITH  CONCENTRATED  LOAD  AT  THE  MIDDLE 

To  show  that  62J  must  be  positive  serai-definite  we  use  the  example  of  a 

simply-supported  beam  with  a  concentrated  load  at  the  middle.  If  the  load  is 

suddenly  removed  (ref  7),  Eq.  (63)  becomes 

0tt  +  b20xxxx  =  0  (81a 

^Shen,  C.  N.  and  Wu,  Julian  J.,  "A  New  Variational  Method  for  Initial  Value 
Problems,  Using  Piecewise  Herraite  Polynomial  Spline  Functions,"  ARO  Report 
81-3,  Proceedings  of  the  1981  Array  Numerical  Analysis  and  Computers 
Conference,  1981. 

^Jacobsen,  Lydkis  and  Ayre,  Robert  S.,  Engineering  Vibrations,  McGraw-Hill, 
1958. 


Or  from  Eq.  (56)  we  have 


ytt  +  h^xxxx  3  0 


(31 


The  starting  conditions  are 


d 

—  o  (x)  =  0 
dt  ° 


oQ(x)  -  Me/ I 


(Px/2)h 


for  0  <  x  <  1/2 


and 


P(£/2)h  x 

on(x)  - - (1  -  ~)  for  1/2  <  x  <  1 

I  £ 


The  solution  for  Eq.  (31b)  Is 


92y 

o(x,t)  -  -  Eh  -~5 
dX 


6a 


where 

and 


>ag  00  1 

--  l  --  (-l)(n-l)/2  3tn(nTTx/£)cos  pnt 

n=odd  n 

pn  =  bn2ir2/£2 


os  =  (P£/2l)(h/2) 

The  quantity  og  is  the  Initial  static  stress  at  the  middle  of  the  beam 
x  *  1/2  and  on  the  outer  surface  of  the  beam. 

In  order  to  find  yt  we  let 

n-1 


1  8os  *  1  A  ,  " 

y  -  -  l  (-oX")  <*!)  2 

n=odd  n  n1T 


-  ^  \  0/\  I  \  * /  sln(nnx/ £)cos  pnt 

Eh  it2  -  •  -z  -« 


Then  by  partial  differentiation  we  have 


a2y  1  80s  * 


n-1 


OOg  ”  1  - 

— - --  l  (--)(-l)  2  sin(nTTx/ £)cos  pnt 

Eh  n*odd  n 

which  agrees  with  Eq.  (83a),  and 


n-1 


3y  1  8os  °°  1  - 

—  - - --  T  (-■?)  (-1)  2  b  sin(nitx/£)  sin 


(82a) 

(82b) 

(82c) 

(83a) 

(33b) 

(83c) 

where 

(34a) 

(34b) 


a.. 


2 


(84 


where  from  Eq.  (83b) 


b  -  pn£2/(n2ir2) 


(84d) 


bir 2/  £2 


(84e) 


It  is  noted  that  the  index  n  appears  in  both  spatial  and  temporal  functions  in 
Eqs.  (84a)  and  (84b)  under  the  summation  sign.  We  are  interested  in  finding 
those  functions  of  t  that  are  Independent  of  the  index  n.  Let  us  assume  that 

n2m 

pnt  =  ---  (A/2-c) 

JC 

*  n2m/2  -  n2irc/&  (85) 

It  is  noted  that  for  n  *  1,  3,  and  5,  n2n/2  becomes  n/2,  4x  +  ir/2,  and  12n  + 
ir/2,  respectively. 


Thus  we  have 


cos  pnt  *  cos  [n2n/2  -  n2nc/£] 

*  cos  {n/2  -  n2mc/ l\ 

*  sin  (n2*/c/&) 


(86a) 


sin  pnt  =*  cos  (n2irc/£) 


(86b) 


Moreover,  for  c/(£/2)  *  1,  1/2,  and  0 


cos  pnt  31  sin[  (ir/2)c/(  1/2)  ]  ■  1,  0.707,  and  0,  respectively  (87a) 

and 

sin  pnt  *  cos[(it/2)c/(  i/2)  ]  *  0,  0.707,  and  1,  respectively  (87b) 
The  above  functions  are  independent  of  index  n  at  those  values  of  c/(£/2). 
Thus,  Eq.  (84)  may  be  rewritten  at  those  values  as 


32y  *  °s 

---  (x,t-t)  - - sin(Ttc/*)y0(x) 

3x2  Eh 


(88a) 


and 


3y  A  °9 

—  (x,t»t)  *  -  coa( nc/ £)y0(x) 

.  TIL 


(88b) 


where 


n-1 

00  - 


y0(x)  "  I  (-1)  2  sin(mtx/ i)  (88c) 

n*odd 

The  series  terms  in  Eq.  (88c)  are  the  result  of  an  expansion  of  a  triangular 
deflection  of  the  form 


yQ(x)  -  x(A/2) 
yQ(x)  -  2  -  x/(A- 2) 


for  0  <  x  <  i/2 
for  i/2  <  x  <  i 


(88d) 

(88e) 


as  shown  in  Figure  2. 

For  the  images  of  Eqs.  (84b)  and  (84c)  the  time  dependent  terras  become 

cos  pnt(T-t)  -  cos(pnT-pnt)  (39a) 

and 

sin  pn(T-t)  -  sin(pnT-pnt)  (89b) 

The  term  pnT  can  be  obtained  from  Eqs.  (84d)  and  (C'e)  as 

pnT  »  (bn2n2/A2)T  -  n2piT  (89c) 

For  computation  by  the  finite  element  method  the  Increment  in  time  is 
taken  as  T  which  is  defined  as 


A 

T  "  fcb  -  t0  *  (m/pl)(*/2)  (9ua) 

where 

m  -  1,2,3...  (90b) 

Then  with  the  aid  of  Eqs.  (85),  (89c),  and  (90a)  we  have 

pn(T-t)  »  mn2(ir/2)  -  [n2n/2  -  n2nc/£] 

-  (m-l)n2(w/2)  +  n2c/ i  (90c) 

Then  for  the  case  when  ra  ■*  l,  the  time  dependent  terms  become 

cos  pn(T-t)  -  cos(n2nc/A)  (91a) 
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and 


(91b) 


sin  pn(T-t)  =»  sln(n2xcM) 

By  similar  method  we  can  obtain 

32y  “  aS 

— r  (x,t*T-t)  - cos(Ttc/ fc)y0(x) 

3x2  Eh 

3y  *  k°s 

—  (x,t»T-t)  - sin(ncM)y0(x) 

3t  Eh 


(92a) 

(92b) 


For  the  partial  differential  equation  we  have  the  second  variation  from  Eq. 
(52)  which  gives 

ip  A  A 

62J  -  /  [<5yt(x,t»t)6yt(x,t=T-t) 

o 

A  A 

+  b26yxx(x,tJ»t)  iSyxx(x,t=“T-t)  ]dt  (93) 

Separating  Eq.  (93)  Into  two  parts  and  using  Eqs.  (88)  and  (92),  we  have 

62J  -  62J[6yt]  +  6ZJ[b6yxx]  (94) 

The  first  terra  on  the  right  of  Eq.  (94)  is 

xb  piT*n/2 

fi2J[Syt]  2  /  (bos/Eh)2y02(x)d xj  cos(  xc/*)sin(  xc/  X.)d(pit) 


where 


■  J* 

xo 

(bog/Eh)2y02(x)dx  >  0 

(95a) 

Pit  -  x/2  - 

xc/*  d(pit)  ■  -(x/*)dc 

(95b) 

Pit  -  x/2 

,  c/(*/2)  ■  0  ,  xc/ *  »  0 

(95c) 

Pit  -  0  , 

c/( Jt/2)  -1  ,  ire/*  -  ir/2 

(95d) 

lq.  (94)  is 

*b 

(bo8/Eh) 

n/2 

2y02(x)dx/  sin(xc/*)  cos(  xc/ *)d(  xc/ *) 

xo 

o 

xb 

-  J 

(bo8/Eh)2y02(x)dx  >  0 

(95e) 
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*/(2pi)  *  t/4 


(99b) 


T  “  tb  -  tc  - 

If  m  “  2  and  S2.l  ■  0  In  Eq.  (99a),  we  can  conclude  that 

T  -  tb  -  t0  <  n/Pl  -  t/2  (99c) 

This  Is  the  upper  limit  of  the  Increment  we  choose  for  T,  above  which  the 
minimum  of  the  functional  J  Is  not  guaranteed. 

CONCLUSIONS 

The  functional  In  bilinear  form  Is  symmetrical  about  the  original 
variables  and  the  adjoint  variables.  The  Euler-Lagrange  equations  for  the 
original  and  the  adjoint  systems  are  derived  using  the  fundamental  lemma  of 
the  calculus  of  variations.  By  integrating  the  bilinear  expression  by  parts, 
one  can  obtain  the  bilinear  concomitant  In  terms  of  Initial  and  boundary 
terms.  The  adjoint  system  can  be  arranged  In  a  manner  that  It  Is  a  reflected 
mirror  of  the  original  system  In  time.  Thus  the  Initial  conditions  for  the 
bilinear  concomitant  become  zero. 

Generalized  boundary  conditions  using  many  types  of  "springs”  relating 
the  various  spatial  partial  derivatives  are  defined  to  satisfy  the  boundaries 
of  the  concomitant.  The  higher  partlals  In  original  variables  and  variations 
in  the  adjoint  variables  can  be  kept  In  low  orders  by  these  "springs". 
Algorithms  are  developed  for  use  in  the  finite  element  method  by  taking  the 
first  variations  of  the  functional.  These  algorithms  are  simplified  because 
the  adjoint  system  gives  exactly  the  same  solutions  as  that  of  the  original 
system. 

Sensitivity  coefficient  Is  found  to  be  similar  to  the  variation  of  the 
variable  and  obeys  the  same  partial  differential  equation.  The  solution  of 
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the  original  PDE  is  taken  as  the  solution  of  the  variations  for  two  examples, 
a  simple  oscillator  and  a  simply-supported  beam  with  load  at  the  middle.  It 
is  found  that  the  second  variation  of  the  functional  is  positive  semi-definite 
if  the  Increment  in  time  for  the  finite  element  method  is  less  than  half  the 
natural  period  of  the  physic  systems  in  both  cases.  This  will  guarantee  a 
minimum  for  the  functional  and  thus  the  method  is  truly  workable  if  employed 
as  algorithms  for  the  finite  element  method. 
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